/*
 * EquilibriumConditions.cpp
 *
 *  Created on: 11 Nov 2011
 *      Author: allan
 */

#include "EquilibriumConditions.h"

EquilibriumConditions::EquilibriumConditions()
{}

EquilibriumConditions::EquilibriumConditions(const vector<EquilibriumCondition>& conditions)
{
	SetConditions(conditions);
}

void EquilibriumConditions::SetConditionMatrix(const MatrixXd& alpha)
{
	this->alpha = alpha;
}

void EquilibriumConditions::SetConditionFunctions(const vector<FunctionTPna>& hj)
{
	hj_function = hj;
}

void EquilibriumConditions::SetConditionFunctions(const VectorXd& hj)
{
	hj_function.clear();
	
	hj_constant = hj;
}

void EquilibriumConditions::SetConditions(const vector<EquilibriumCondition>& conditions)
{
	// The number of equilibrium conditions
	const uint Nj = conditions.size();
	
	// The number of species handled in the equilibrium conditions
	const uint N = conditions.front().alpha_row.rows();
	
	// Allocate memory for the matrix alpha and container hj_function
	alpha.resize(Nj, N); hj_function.resize(Nj);
	
	// Loop over all the equilibrium conditions
	for(unsigned j = 0; j < Nj; ++j)
	{
		// Set the j-th row of the alpha matrix
		alpha.row(j) = conditions[j].alpha_row;
		
		// Set the j-th row of the equilibrium condition functions
		hj_function[j] = conditions[j].hj;
	}
}

const MatrixXd& EquilibriumConditions::GetConditionMatrix() const
{
	return alpha;
}

const VectorXd EquilibriumConditions::EvaluateConditionFunctions(double T, double P, const VectorXd& n, const VectorXd& a) const
{
	if(hj_function.empty())
		return hj_constant;
	else
	{
		const uint Nj = alpha.rows();
		
		VectorXd hj = VectorXd::Zero(Nj);
		
		for(uint j = 0; j < Nj; ++j)
			hj[j] = hj_function[j](T, P, n, a);
		
		return hj;
	}
}
